BROWN-HET-1360; CU-TP-1089; MIT-CTP-3389 

June 2003 
hep-lat/0306038 



A Review Of Two Novel Numerical Methods in QFT 



R. Easther 1 , D. D. Ferrante 2 , G. S. Guralnik 3,4 and D. Petrov 2 



1 Department of Physics 
Columbia University, NY — 10027, NY. 



2 Department of Physics 
Brown University, Providence — 02912, RI. 

3 Center for Theoretical Physics 
Massachusetts Institute of Technology, Cambridge — 02139, MA. 



Abstract 

We outline two alternative schemes to perform numerical calculations in quantum field 
theory. In principle, both of these approaches are better suited to study phase structure 
than conventional Monte Carlo. The first method, Source Galerkin, is based on a numerical 
analysis of the Schwinger-Dyson equations using modern computer techniques. The nature 
of this approach makes dealing with fermions relatively straightforward, particularly since we 
can work on the continuum. Its ultimate success in non-trivial dimensions will depend on the 
power of a propagator expansion scheme which also greatly simplifies numerical calculation 
of traditional perturbation graphs. The second method extends Monte Carlo approaches by 
introducing a procedure to deal with rapidly oscillating integrals. 
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1 Introduction 



Over the last two decades, Monte Carlo numerical methods applied to lattice quantum field theory 
have allowed us to perform many important calculations, and have verified and extended our basic 
understanding of elementary particle phenomena. 1 While very informative, these calculations still 
do not fully reflect the potential of these methods because it has been impossible to obtain the 
multiple teraflop- years of computer time required for very accurate non-quenched fermionic calcu- 
lations. The promise of "definitive" lattice QCD calculations has helped to drive the development 
of supercomputers to the point where computer clusters are beginning to sustain speeds of several 
teraflops. Over the next few years, many teraflop-years will be devoted to lattice QCD, and it is 
likely that the range and accuracy of prediction of QFT using numerical calculation will be greatly 
extended. 

However, even with this promise of great resources, numerical field theory calculations are far 
from becoming mechanical. While large blocks of computing time on fast new machines should 
make it possible to deal with the fermion determinant, the effort of doing so and dealing with lattice 
artifacts guarantee that it will still be difficult to extract reasonable numerical information about 
fermionic systems. Moreover, calculations involving non-positive definite actions, actions with 
rapid oscillations, phase transitions and symmetry breaking are generally resistant or intractable 
to solution by Monte Carlo methods. Many problems, such as general scattering, are just too 
computationally intensive to be accessible. Motivated by these issues, we are in the process of 
developing two supplementary approaches to traditional numerical quantum field theory. The first 
is the "Source Galerkin Method" (SG) 013111 which uses nested approximations to the Schwinger- 
Dyson equations. The second is a tuned Monte Carlo method, "Mollified Monte Carlo" (MMC) [HI 
Hm E] which uses an averaging process to smooth regions of high oscillation in path integrals. 
Conceptually distinct, these methods actually have similar starting points. MMC uses information 
from stationary phase points of an action, while SG methods are easiest to implement if they are 
iteratively constructed around stationary phase points. 

2 Source-Galerkin Methods 

The Monte Carlo approach suffers from at least two serious intrinsic problems associated with 
fermions. The first is that, in order to apply this technique, it is essential to formulate a field theory 
on a space time lattice rather than on the continuum. Any calculation, of necessity, involves an 
extrapolation to the continuum for the final result. As a consequence of the lattice, it is necessary to 
introduce extra degrees of freedoms for fermionic fields. This produces artifacts which are difficult 
to control. The second problem arises from the fact that we can not directly perform Monte 
Carlo sampling on integrals with the anti-commuting Grassman variables requisite for fermions. 
Consequently, any fermionic degrees of freedom must be explicitly removed by partial integration 
of the lattice path integral. These integrations result in the very non-local "fermion determinant" 
which requires multi-teraflop computer power to evaluate to high accuracy. 

We introduced the Source Galerkin method to deal with these two fermionic problems and 
the difficulties associated with phase structure mentioned in the introduction. The SG method is 
based on an iterative set of approximations to the Schwinger-Dyson equations. It is defined on the 
continuum and has the nice property that fermions (except for anti-commutativity) are treated 
in the same way as bosons. Most of the basic concepts of the SG method can be demonstrated 
by its application to single field models. In particular, we can get some understanding of why it 

This document is based on a talk given by G. Guralnik at the "Seventh Workshop on Quantum Chromodynam- 
ics", 6-10 January 2003. 
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is possible to deal with phase structure issues. Therefore, in the interest of simplicity, we avoid 
directly addressing the details of fermionic calculations here and confine this discussion to these 
models. The basic ideas of how fermions work can easily be extrapolated from a primitive lattice 
precursor of this model [5]. 

We start with the simple assumption that the Quantum Field Theory (QFT) action is written 
with sources Ji(x) for every field 4>i(x) so that the vacuum functional Z = (+010—)^ satisfies the 
differential equation, 

essentially the 
field equations 

A familiar but far from straightforward example of this is 4 scalar QFT 

(d 2 + m 2 ) (f)(x) + g (f) 3 (x) = J(x) 



which, in Euclidean space, gives the very non-trivial set of coupled differential equations: 



{d 2 x + m 2 



5J(x) * \5J(x) 



5 \ 3 



Z[J] = 



2.1 The Ultra-Local <\y 




To further simplify, we examine the above quartic field theory in 
zero dimensions. We will show that even this case has a rich basic 
structure which is not directly accessible by normal Monte Carlo 
approaches. Our demonstration example then reduces to the linear 
differential equation 



m "dJ +9 \dJJ 



Z[J] = . 



(1) 



Figure 1: Path integrals 

must start and end at oo in Since this has three independent solutions, three independent pa- 
rameters are required to specify the full solution pQ. If we rewrite 
Z[J] as a power series, 



hatched regions. 



oo 1 

Z[J} = Z[0}Y,^G k J k ; 



k=0 



equation ((TJ leads to the recursion relations: 

g G n+3 + m 2 G n+1 - n G n _i = 



These have three independent solutions that can be expressed in terms of parabolic cylinder 
functions. pQ This is a useful test system, since it is computationally non-trivial, but has solutions 
in terms of well known functions, making it is easy to check the validity of a new numerical 
approximation technique. Note that, with the proper choice of parameters, solutions exist with 
non-vanishing odd Green's functions. 



3 



Positive 




Imaginary 






Real 


— oo 


+00 


Negative 




Imaginary \ 


^ —zoo 



Figure 2: Three convenient 
independent paths. 



It is usual to think of the path integral formulation of QFT as 
being defined by a process involving evaluation of integrals along 
the entire real axis. However, this cannot yield a complete de- 
scription of the ultra local quartic theory, since it is impossible to 
generate non-zero odd-order Green's functions from actions with 
only even powers of fields, when path integrals are only evaluated 
in this way. As a generalization of this, we see that a choice of 
path confined to the real axis removes the possibility of sponta- 
neous symmetry breaking! Very specifically, it is straightforward 
to show that the one dimensional path integral along the whole real 
axis for the above zero dimensional field theory does not produce 
the full solution set outlined above. Fortunately, it is easy to con- 
struct the path integral formulation so that it coincides with the 
Schwinger-Dyson approach 0. For our simple model, the "path 
integral" 



Jf TC? 
expj- — 



x 



9 4 



JX 



dx 



satisfies the Schwinger-Dyson equations as long as the integrand vanishes (or the exponent goes to 
—00) at the ends of the integration region. This happens if the integration begins and ends inside 
the cross hatched regions shown in figure 1. It follows, in agreement with our differential analysis, 
that three independent paths can be chosen in this model, yielding three independent solutions. 
A convenient choice of paths is: 

• Real: f G E; 



Positive Imaginary: £ G 



Negative Imaginary: £ G 



-oo,0] U [0,+zoo 
-oo,0] U [0 



I 00 



These choices are displayed in figure 2. 

Note that, because we are solving a third order linear differential equation, the "zero-field", "one- 
field" and "two-field" Green's functions are arbitrarily determined by the choice of boundary condi- 
tions and not by the dynamics! It is convenient to group solutions together in 3 real combinations 
which can be characterized by their degree of singularity around vanishing coupling: 



• Regular at g — > 0: consistent with perturbation theory; 

• ~ 9~ 1 ^ 2 at g — > 0: "symmetry breaking"; 

• ~ exp(/z 2 /4g) at g — > 0: "instanton". 



It can easily be shown that the behavior of these solutions around zero coupling are related to 
expansions around the three stationary phase points of the path integral in the complex £-plane. 
Indeed, direct calculation shows that all three of these expansions form asymptotic series corre- 
sponding to expansions of the exact solutions. 
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2.2 Implementing the Source- Galerkin Approach 

An approximate solution to the differential equation 

D[u] = , (2) 
can be constructed as a linear combination of the first N members of a set of trial functions {</)}. 

TV 

Now we can develop a procedure, which will enable us to determine a set of coefficients {a} so that 
the solution u a defined by them is as "close as possible" to the exact solution for a particular choice 
of {0} and N. We start by obtaining the residue R of an approximate solution by substituting the 
u a into the equation (J2J). 

R{a,x)=D[u a ]^0. 

We define a scalar product of two functions and choose a set of weight functions {w}. Next, a set 
of equations for determining ai is constructed by requiring that the scalar product of the residue 
with the weight functions vanishes: 

(R\wi) = 0, Vi = 0,...,JV. 

This approach is an example of "The Method of Weighted Residues". Provided that the Wi form a 
complete set of functions, it guarantees that lim/v^ooW a = u in the mean. If the specific choice of 
Wi = (pi is made, then this approach is referred to as "The Galerkin Method". 

There are several points that need to be taken in to account when using this method to solve 
field theories. As stated above, the Galerkin technique produces approximate solutions which 
converge to the exact solution as the number of terms goes to infinity. In practice, the resulting 
series has to be truncated. Our experience to date suggests that for nontrivial theories, calculations 
for N > 8 will be impossible to perform. Consequently, when we pick functions to describe our 
problem, it is essential to make a reasonable guess consistent with all the expected properties and 
symmetries in order to achieve rapid convergence, even with few terms. In practice the apparent 
single pole dominance of many important qualities in physical theories, which is so important 
to Monte Carlo success, seems to serve us equally well in the Source Galerkin approach. Finally, 
picking a definition of an inner product which is efficient to use with the chosen expansion functions 
greatly facilitates rapidly accurate computations. 

2.3 The Actual Numerical Solution 

In order to test the validity of the SG calculational technique, we applied it to the ultra-local 4 . 
We solved this model using two different sets of trial functions: powers of J and a set of Hermite 
polynomials in J. The approximate solutions constructed from these trial functions are of the 
form: 

N . N 

z a [j] = Y.j\ akJk ' Z »V] = Y1 ak H ^ ■ 

k=0 ' k=0 

The solution in terms of truncated power series is easy to deal with and maps directly to the 
theoretical solution obtained in section 12.11 However, these are not orthogonal functions and 
simple generalizations for more complicated problems do have some difficulties. On the other 
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hand, Hermite polynomials are orthogonal under proper choice of inner product. This property 
greatly simplifies the process of determination of the coefficients a, and ensures high numerical 
stability 

The residual equations for the truncated power series were obtained using one of the inner 
products presented below: 

1. Canonical: 

€ 

(f\g) = j f{J)g{J) dJ; 

—e 

2. Exponentially weighted: 

oo 

(f\g) = J f(J)g(J)e- j2 ^ 2 dJ; 

— oo 

Note that both of these inner products are chosen in such a way that they emphasize a small 
region close to zero. This is done because the generating function and its derivatives will be 
used to compute the propagators which are evaluated with the sources set to zero. Epsilon is a 
parameter, which allows us to tune the method to test it's stability and achieve better numerical 
accuracy. It becomes highly important in more interesting models. 

Figure El shows dependence of absolute error in the determination of the generating functional 
relative to the source j. It is clear that the accuracy of one part in 10 7 is achieved and exceeded 
as the value of the source approaches zero. 




Figure 3: Absolute error in determination of Z using a truncated power series. 

When using Hermite polynomials as a set of trial functions, the most sensible definition of inner 
product is the one under which the Hermite polynomials are orthogonal to each other. 

oo 

(f\g)= j f(J)g(J)e- j2 ^ 2 dJ. 

— oo 

The total accuracy of the computation given by this choice of trial functions is comparable 
to the accuracy of the solution in terms of the power series and is limited only by the numerical 
precision of the calculation. Order by order the relative error in determining coefficients Oj with 
respect to epsilon is shown on figure 0] Both graphs show the same solution at different scales. 
Clearly the precision of the computation of ai decreases at higher orders. From the graph on the 
left it is evident that a 7 is determined with an error of 1 per cent while the graph on the right 
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shows that a 5 can be calculated more accurately than 5 parts in 10 4 . This trend does not present a 
problem since the absolute value of the contributions decreases at high orders thus compensating 
for the loss of accuracy. 

Note that the solution is constant in a significant range of values of the parameter. At small 
values of epsilon this range is limited by numerical precision of the hardware used for the compu- 
tation. From the graphs it is clear that the calculation becomes numerically unstable when epsilon 
is set to any value smaller than 0.24. For large values of the parameter, the range of stability 
is limited by the number of terms in the approximate solution. If epsilon is set too high, the 
generating functional can not be accurately approximated at every point in the range (— e... + e) 
by the linear combination of some number (seven in this example) of Hermite polynomials. 

The stability property described above is important in computations where a theoretical so- 
lution is not available. Then we tune the parameters of a problem to find a stable region which 
should indicate that we have found a good solution. 




a s 



Figure 4: Relative error of Hermite polynomial expansions. 



Now, we return to higher dimensions. In order to solve, 

through numerically approximation, we start with the exact answer in the form: 

{/ G 2 J 1 
a + J(x) G(x) + — -— + J J J G 3 + J J G 4 J J + • ■ ■ \ . 

Here, co-ordinates are only displayed in the linear term and space-time integrations are not shown. 
Initially, the G l ( be anything consistent with the symmetries of the theory. The 

fact that the zero dimensional model did not have integrations is why it was so simple. 

Using the Poincare invariance to produce spectral forms, we conjecture that it should be possible 
to span the solution space (or a very large part of it) with products of arbitrary two field propagators 

G P (x - y) = f ^^^^f^ d% d 2 K . 

Thus, we represent G n by all possible graphs with n external lines, all numbers of lines intersecting 
at any vertex and a p (K 2 ) different for each propagator. Important! This is NOT PERTURBA- 
TION THEORY. The masses and weights in each of these propagators are arbitrary at this stage. 
If our conjecture is correct, we have merely exploited the restrictions of relativity and translation 
invariance in writing this representation. The structure we have written is still extremely complex. 

The above expansion inserted into the field equations yields conditions of constraint on a p (K 2 ). 
In general, they are too complicated to solve directly. After all, field theory has not be exactly 
solved! To proceed, we simplify the problem as follows, 
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• Truncate the expansion in J; 



• Limit the number of masses in each propagator; 

• Limit the number of graphs considered for each J. 

We do this in an organized, systematic way, so that after a first guess more terms can be included 
to iterate on the initial answer. Roughly, our procedure works as follows: Begin by guessing an 
initial solution Z approx [J] consistent with space-time, spin and other internal symmetries. Since we 
have not guessed an exact solution: 



The idea of Source Galerkin is to require that 



J dJ\ ■ • ■ dJ n i^j(<7i, . . . , J n ) F ^^^j^j ^approx 



[J]=0 



so that % a p prox [J] satisfies the field equations "on the average". The number of (in principle arbitrary 
functions) Fj is picked so that all the undetermined weights in 2, a p prox [J] are determined. 

In general, the equations to be solved are non-linear. Consequently, solutions must be de- 
termined in a very careful and systematic manner and parameters must be tuned. Theorems for 
simpler Galerkin-like approaches promise convergence and we assume this is true for this approach. 
We have tried this method in many simple models, with amazing accuracy when a check is avail- 
able. In general, it must be confirmed that the results are stable and convergent. Higher order 
iterations, while simple in principle, have caused us computational difficulties in practice. While we 
believe that these problems can be resolved, we have not yet fully demonstrated that our method 
is practicle in more interesting theories. 



2.4 The Nonlinear Sigma Model 

Here we give a simple but interesting example of a lowest order SG solution. Calculation to the 
accuracy obtained in this case is moderately difficult to replicate in Monte Carlo. We have also 
examined leading order corrections which, as stated, are essential for the proof of validity of the 
method, and which we will discuss elsewhere. The theory is given by 



l - I d, &> + X {X) (W) a (x) - p 



After the introduction of sources J a (x) and S(x) for canonical and auxiliary fields, one has 



1 

V 9 



£ = ^3„ <j> a {x) V a (x) - J a (x) a (x) + x {x) ( a (x) <f> a {x) - - ) - x(x) S(x) 



The equations of motion that follow are: 



- \z[J,S] -2S{x)Z[J,S] = . 



5J a (x)SJ a (x) g 2 
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Z a [J,S] = Z expji J 



We guess the leading approximation to be 

The residues are computed by substituting this expression into the Schwinger Dyson equations, 

Rx = □ J G a \x - J\0 d£ - J a (x) -xoj G a \x - J 6 (£) di 

R 2 = G aa (x - x) + J G a \x - J b (0 d£ J G a \x - rj) J b (7]) d7]-\-2 S(x) . 

The Source Galerkin equations can be obtained by requiring that the projections of the residues 
on the weight functions vanish: 

<l|<Ay)| J R 1 (x)) J > 5 = 0; 
(11(11^)^ = 0; 

which yield, 



□ G ab (x -y)- 5 ab (x -y)- X o G ab (x - y) = ; 
G aa (0) + e 2 J G ab (£ - V ) G ba ( V - £) - \ = ; 

where, the following definition was used for the inner product 



(A\B)j = ±-J A[J(x)]B[J(x)]e- j2 ^ 2 [dJ] 



It is straightforward to show that the results obtained here give the same two field Green's function 
as the leading order large N expansion. 

As indicated, the problems with our method occur at higher orders in the expansion. Since 
our approach can be represented by arbitrary combinations of Feynman graphs whose masses and 
coefficients are set by solving Galerkin equations, fast numerical evaluation of these graphs is 
essential for Source Galerkin to be useful. It is for this reason that we developed a new method to 
numerically study Feynman graphs. 

2.5 Feynman Diagrams 

Motivated by our work on the Source Galerkin technique, we have carefully looked at the numerical 
calculation of Feynman graphs. Historically complicated graphs are evaluated with Monte Carlo 
methods which take a relatively long time to achieve accurate results. We have have been able 
to devise a new method for numerically evaluating graphs which is very accurate and very fast, 
and has considerable promise for complex traditional perturbation calculations (such as the eighth 
order magnetic moment). The availability of such a method is crucial to fully implementing our 
Galerkin approach. 

We use an approximation to the propagator that reduces graph evaluation to a rapidly conver- 
gent mult i- dimensional sum.jH] We construct a propagator representation as follows: The cutoff 
two field propagator, G\ 

(2 7r) 4 p 2 + m 2 
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is written using the Sine function, 



as 



m 2 h 



G A , h {x) = j—^ Yl ex P 



fc= — oo 



2 2 

m x 



AC(k)_ 



_2 fcfc-e" 
kh m ' ^ e 



A 2 C 2 (A;) 

The accuracy of the approximation is tuned by the value of /i, and typically, the propagator can 
be approximated to very high accuracy (1 part in 10 16 ) with fewer than 100 terms in the sum. 
Inserting this representation of the propagator into the graphical topologies we wish to calcu- 
late ensures that all the space-time or momentum integrals are Gaussians and can be performed 
with ease. The resulting multi-dimensional sum can then be quickly and accurately computed 
numerically. We have calculated 4th non-trivial order (3 loops) in scalar field theory [H] and suc- 
cessfully responded to challenges to the accuracy and speed of the method. We have computed 
the QED graphs shown in the diagram on the right [7] as well as a few higher order graphs. 



Technical issues with accuracy using an "auto" renormalization 
scheme have slowed the development of a fully automated imple- 
mentation of this method to higher order QED graphs. We believe 
these difficulties are soluble. Moreover, the power of this graphical 
calculation method gives us reason to believe that we can iterate 
the SG method to moderate order with relatively small amounts of 
computer time. If this is true, a much broader range of QFT will 
be opened to numerical solution with current computer power. 

Finally, we note that there is a long history of trying to solve QFT by making approximations 
to the Schwinger Dyson equations. Our approach is unique because of our way of picking best fits 
and our systematic procedure based on the powerful way of evaluating arbitrary graphs using the 
Sync expansion techniques described above. The current power of computers has made it possible 
to revisit these ideas and institute algorithms which previously would have been useless. 




3 Mollified Monte Carlo 

The success of Monte Carlo methods is very dependent on the form that the "action" takes. If 
some terms are not positive (sign problem) or imaginary, or if the action has regions of very rapid 
oscillation, conventional Monte Carlo approaches will mostly fail. The failed actions are amongst 
the more interesting ones that we would like to evaluate. As a concrete example, consider symmetry 
breaking such as was discussed in the SG sections. We saw that to get non-vanishing odd Green's 
functions, it was necessary to extend the path integral onto the imaginary axis or, equivalently, 
to examine actions with negative bare mass terms and half line integrations. In general we would 
not expect Monte Carlo to be able to handle this. Further it would be very nice if we could, at 
least in some formal manner, examine actions in Minkowski space, which means we must deal with 
(carefully defined) complex integrals. Again, we can not do this conventionally. The study of phase 
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structure means looking at actions in regions of high oscillation. Using normal approaches, this is 
asking for trouble. 

Doll and his collaborators jHl Ell have discovered a way of rewriting oscillatory integrals which 
has the potential to overcome most of these problems and allow Monte Carlo like numerical eval- 
uation. We will outline some of the fairly complex details of this approach in what follows, but 
it is straightforward to state the general procedure. Basically, we will show by the use of a prob- 
ability function it is possible to rewrite (mollify) an oscillating path integral, so regions that do 
not contribute (regions of non-stationary phase) are effectively removed. Next we sample the new 
integral form using importance functions which are heavily weighted around the stationary phase 
regions of the mollified integral. This modified Monte Carlo approach will converge even with in- 
tegrals which originally possessed awful oscillations. In fact, Mollified Monte Carlo (MMC) tends 
to yield results near to mean field (large N) analytic approaches. The results will also be close to 
SG expansions picked in our usual way. We have examined this "mollified" approach extensively 
in zero dimensional field theory and are now trying to use it to solve real theories. While we see 
rapid convergence in regions of high oscillation, this method is more complex than normal Monte 
Carlo and will not save any time when used in non-oscillatory problems. Further, it provides no 
new insight into unquenched fermionic calculations - only the potential of convergence in regions 
previously not accessible. 

We make this more specific by showing some simple examples. Our starting point is a seemingly 
innocent identity using a probability distribution, We define: 



f(x)dx» / (f(x)) e dx, (3) 



where 



(f(x)) e = (f*P t )(x), (4) 
Pe(x-y)f(y) dy, 



with 



P e {x) dx = 1 



In this case, the smoothing function is given by P e (x) which is called a "mollifier" or an "approximate 
identity". We use one of the properties of the mollifiers, 



(f(x))/-^f(x) / (f( x )) e dx^ I f(x)dx. 



With this simple averaging method we can build mollified versions of the desired partition functions 
and Green's functions, tame oscillations and evaluate previously numerically inaccessible functions. 
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3.1 Simple example 



Let us show how oscillations can be tamed by mollifying the very singular function exp(i/x) 
cos(l/x) +i sin(l/sc). 



g(x) = exp(i/x) , P e {x) 



exp 



2 e 2 J 



v/2 



7re z 



exp{-(l/2) (x - y) 2 /e 2 } exp{i/7/} 



V2 7T6 2 



0.2 0.4 0.6 0.8 1 



Re[mollg(x,e=1)] 



0.2 0.4 0.6 0.8 1 



An analogous behavior is seen for the imaginary part. It is quite clear that we have smoothed an 
essential singularity We will show that this can be used without loosing information to evaluate 
Green's functions. 



3.2 Simple example (3D): Gaussian mollifier 

We give one more example to show that the above result was not a unique case: 



f m (x) = exp (i ma; 2 ) , P e (x) 



-If 



V2 



(fm(x)) t 
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Re[f] (red) and Re[mollf] (blue) 



lm[f] (red) and lm[mollf] (blue) 






Once again, it is easy to see that the surfaces in blue/dark grey (mollified) are less oscillatory than 
their counterparts (red/light grey). This is quite useful when one starts to think about combining 
the mollifier technique with the Monte Carlo one and sees that we have tamed this highly oscillatory 
function for positive and negative values of mass and for real and imaginary contributions. This 
makes it clear that we should be able to handle general sign problems. 

3.3 Application of the Method 

Now we can outline the application to Monte Carlo integrals. The method, consists of two basic 
steps. We start by mollifying (pre-averaging) the integrand. The generating functional is given by, 



using the following [Gaussian] mollifier: 

Pe(y) 



e iS(x)-iJ-x 



X\ 



exp<-iy T - (e 2 )- 1 -y 



y/(2ir) n det(e 2 ) 



We next perform a stationary phase expansion which produces the structure we will use for the 
importance sampling functions. We start with the basic integral 



J </(x)e i5(x) " ij - x ) e [dx] 



and derive the result, 



/(x )e i5(xo) - ijxo (5) 
exp{-± B T ■ (1 + e T ■ S" • e)" 1 • b} 



det(l + e T -S"-e) 



[dx] , 

where e 2 is the covariance matrix that defines the Gaussian distribution, B = e • V5(x) and 
[S"]ij = d 2 S(x.)/dxidxj. Note, we have suppressed the important point that, in general, the 
functions under study will have multiple stationary phase points. 
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The integrals of interest take the form: 



r 

I 


/ e i5(x)-ij-x\ 

w(x) 


[dW(x)} 


J 


f 




[dW(x)} 


W(x) 



(6) 



where W(x) is known as the importance function and the MC (Metropolis) sampling is done over 
this new measure. A reasonable properly biased choice is given by, the result of our saddle-point 
expansion. Thus: 



WJx] 



expjz S(x ) - \ B T • (1 + e T • S" ■ e)" 1 ■ b} 



^/dfitp + e T • S" • e) 



3.4 Airy Function 

The action for this model is given by: S(x) = x 3 /3 + tx. This is an interesting model because 
one can explicitly calculate the partition function and compare it with the results coming from 
our mollified Monte Carlo procedure. Moreover, this model has two stationary phase points and 
the results we display are from the one in the complex plane which is not accessible with normal 
Monte Carlo. The quantity of interest is: 



Z[t] 



Ai(0) 



The graphs below show the behavior of the partition function. In the leftmost one, the oscil- 
latory behavior can be clearly seen and, in the rightmost one, the same behavior is shown (for a 
fixed t = —16, in blue/dark grey) together with the mollified version of the partition function (in 
red). 



Plot of Fte(e*(-S(x,t))) 




Kit lii 1 ' 



III! Illl 



It is easy to perform the detailed calculations we outline above to get the pure numerical 
MMC evaluation of this integral. The results are sensationally accurate and show that in this 
simple example that we can achieve results not allowable in ordinary Monte Carlo approaches. 
Further, we have done these calculations for the quartic ultra local theory and achieve similar 
highly accurate results corresponding to all three solutions. 
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Conclusions 



We have examined two possible numerical methods that could extend the results of normal Monte 
Carlo approaches. We have demonstrated that, at least in zero dimensional field theory, these 
methods make sense and have unique and useful properties. We only have glimpses of the validity 
of either method in higher dimensions. However, the development of the higher dimensional SG 
analysis has led to a new and very rapid method of doing normal perturbation theory. We expect 
to be able to push both of our methods considerably further in the next year. However, full 
confirmation of their usefulness will probably have to wait for work by other groups since, like 
all numerical quantum field theoretic calculations, large collaborations will be required for hard 
calculations. 
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